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0) Introduction. 

• We study in this article the problem of boundary conditions for the 
equations of gas dynamics. We restrict ourselves to the model of perfect fluid 
and essentially to simple unidimensional geometry. The problem can be consid- 
ered from two points of view. First the linearization of the Euler equations of 
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gas dynamics allows the development of the so-called method of characteristics 
and the linearized problem is mathematically well posed when boundary data 
are associated with the characteristics that comes inside the domain of study. 
Second the physical approach separates clearly fluid boundary conditions where 
incomplete physical data are given, and rigid walls where the slip condition de- 
scribes the interaction between the fluid and a given object located inside or 
around the fluid. There is at our knowledge no simple way to connect these 
mathematical and physical approaches for strongly nonlinear interactions. 

• We introduce in this contribution the notion of partial Riemann pro- 

blem. Recall that the Riemann problem describes a shock tube interaction 
between two given states ; the partial Riemann problem is a generalization of 
the previous concept and introduces the notion of boundary manifold. In what 
follows, we first recall very classical notions concerning gas dynamics and the as- 
sociated Riemann problem. In a second part, we introduce the partial Riemann 
problem for general systems of conservation laws and proves that this problem 
admits a solution in some class of appropriate nonlinear waves. In section 3, we 
recall the linearized analysis with the method of characteristics, introduce the 
weak formulation of the Dirichlet boundary condition for nonlinear situations in 
terms of the partial Riemann problem and show that lot of physically relevant 
situations are described with this theoretical framework. In the last paragraph, 
we propose a practical implementation of the previous considerations with the 
finite volume method. 



1) EULER EQUATIONS OF GAS DYNAMICS. 
1.1 Thermodynamics. 

• We study a perfect gas submitted to a motion with variable velocity in 

space and time. We note first that the primitive unknowns of this problem are 
the scalar fields that characterize the thermodynamics of the gas, i.e. density p, 
internal energy e, temperature T, and pressure p (see e.g. the book of Callen 
[Ca85]). In what follows, we suppose that the gas is a polytropic perfect gas ; 
it has constant specific heats at constant volume C v and at constant pressure 
C p . These two quantities do not depend on any thermodynamic variable like 
temperature or pressure ; we denote by 7 their ratio : 

(1.1.1) 7 = ^ (= Cste). 

We suppose that the gas satisfies the law of perfect gas that can be written with 
the following form : 

(1.1.2) p = ( 7 -l)pe. 
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As usual, internal energy and temperature are linked together by the Joule- 
Thomson relation : 

(1.1.3) e = C V T. 

• In the formalism proposed by Euler during the 18° century, the motion 
is described with the help of an unknown vector field u which is a function of 
space x and time t : 

(1.1.4) u = u(x, t) . 

In the following, we will suppose that space x has only one dimemsion (x G 
IR). We have four unknown functions (density, velocity, pressure and internal 
energy) linked together by the state law (1.1.2) In consequence, we need three 
complementary equations in order to define a unique solution of the problem. 

• The general laws of Physics assume that mass, momentum and total 
energy are conserved quantities, at least in the context of classical physics asso- 
ciated to the paradigm of invariance for the Galileo group of space-time transfor- 
mations (see e.g. Landau and Lifchitz [LL54]). When we write the conservation 
of mass, momentum and energy inside an infinitesimal volume dx advected 
with celerity u(x, £), which is exactly the mean velocity of particules that com- 
pose the gas, it is classical [LL54] to write the fundamental conservation laws of 
Physics with the help of divergence operators : 

d.1.5) | + £(H=0 

d c) 
(1.1.6) -{pu) + — (pu 2 +p) = 

(LL7) dt(l pu2 + pe ) + ^((^ w2+ ^) w + ^ w ) = °- 

We introduce the specific total energy E by unity of volume 

(1.1.8) E = \u 2 + e 

and the vector W composed by the "conservative variables" or more precisely 
by the "conserved variables" : 

(1.1.9) W = (p, pu, pE) 1 . 

The conservation laws (1.1.5)-(1.1.7) take the following general form of a so- 
called system of conservation laws : 
OW d 

( L1 - 10) Ht + Y X F{W) = 

where the flux vector W i — > F(W) satisfies the following algebraic expression : 
(1.1.11) F(W) = (pu, pu 2 + p, puE + pu) 1 . 

• The system of conservation laws (1.1.10) is an hyperbolic system of 
equations sometimes stated as a "quasilinear system" of conservation laws. From 
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a mathematical point of view, the actual state of knowledge (see e.g. Serre 
[Se96]) does not give any answer to the question of the existence of a solution 
(x, t) i — > W(x, t) when time t is no longer small enough, even if the initial 
condition x i — > W(x, 0) is a regular function of the space variable x. Even if 
the problem is solved in the scalar case (Kruzkov [Kv70]), nonlinear waves are 
present in system (1.1.10), they can create discontinuities, and this phenomenon 
makes the general mathematical study of system (1.1.10) very tricky. Moreover, 
the unicity of irregular solutions of this kind of system is in general not satisfied. 
The usual tentative to achieve uniqueness is associated with the incorporation of 
the second principle of Thermodynamics. Carnot's principle of thermodynamics 
can be explicitly introduced with the specific entropy s. Recall that for a perfect 
polytropic gas, the specific entropy is a function of the state thermodynamic 
variables p and p : 

(1.1.12) s = C„Log(^4) 

and it is easy to see that the following so-called mathematical entropy r/(») 

(1.1.13) rj(W) = -ps 

is a convex function of the conservative variables if temperature is positive (we 
refer e.g. to our proof in [Du90]). It admits an entropy flux £(W) in the sense 
of the mathematical theory proposed by Friedrichs and Lax in 1971 [FL71] and 
we have 

(1.1.14) £(W) = -pus. 

The second principle of increasing of physical entropy when time is increasing 
has been formalized mathematically by Germain-Bader [GB53], Oleinik [0157], 
Godunov [Go61] and Friedrichs-Lax in 1971 [FL71] among others. It takes inside 
the formalism of mechanics of continuous media the following weak "conservative 
form" : 

(1.1.15) m r,{W) + fot(W) < 0. 

Inequality (1.1.15) is exactly an equality for a regular solution W(», •) of con- 
servation law (1.1.10). It has to be considered in the sense of distributions for a 
weak solution W(», •) of conservation law (1.1.10). 

• The mathematical entropy (1.1.13) is a strictly convex function of the 

conservative variables if temperature is positive. We introduce the Frechet 
derivative dr](W) of entropy ?](•) with respect to the conservative variables 
W : 

(1.1.16) drj(W) = (j) dW 

and define the so-called entropy variables (p. The precise calculus of entropy 
variables is elementary from the traditional expression of the second principle 
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of thermodynamics. If a container has volume V, mass M, internal energy S, 
then entropy S is a function £(•) of the above variables : 

(1.1.17) S = S(M, V, 5), 

and moreover, function £(•) is an extensive function (homogeneous of de- 
gree 1) of the preceeding variables, id est E (AM, A V, A£) = XT, (M, V, 6) , 
V A > . The classical differential relation between these quantities can be 
stated as 

(1.1.18) d£ = TdS - pdV + pdM . 

Notice that the third intensive thermodynamic variable \i is just the specific 
chemical potential, i.e. the chemical potential by unit of mass. Then the 
global values are related for corresponding specific quantities with the help of 
the following relations : 

(1.1.19) M = pV, 8 = eM, S = sM , 8 + pV-TS = pM . 
Due to the homogeneity of function £(•), the mathematical entropy ?](•) can 
be simply expressed with the help of this thermostatic primitive function : 

(1.1.20) " ri(W) = -E(p, l,pe). 

Taking into account the relations (1.1.8) and (1.1.19), we easily differentiate the 
relation (1.1.20) and it comes 

1/ \ 7/ I 

(1.1.21) d V (W) = -(p--)dp + -d(pu) - -d(pE). 
By comparison between (1.1.16) and (1.1.21) we have 

(1.1.22) = -(/i- - lU ,-l). 
1.2 Linear and nonlinear waves. 

• In this section, we construct particular solutions of the so-called Riemann 
problem. First recall that the Riemann problem consists in searching an entropy 
solution (i.e. a solution satisfying (1.1.15) in the sense of distributions, see e.g. 
Godlewski-Raviart [GR96] or Serre [Se96]) W(x, t) (x e H, t > 0) of the 
following Cauchy problem : 
dW d 

^ -w + = 

(1.2.2) W( X ,0) = {^ HI 

We first remark that conservation law (1.2.1) is invariant under a change 
of the space-time scale. Space-time transform T\ parameterized by A > 0, 
and defined by the relation T\(x, t) = (Xx, At) can be applied on any (weak) 
solution W(», •) of system (1.2.1)-(1.2.2) and generates a new weak solution 
T x W{x, t) defined by the condition (T A W){x, t) = W{T x {x, t)) . We remark 
also that the initial condition (1.2.2) is invariant under space dilatation, i.e. 



6 



Partial Riemann problem, boundary conditions, and gas dynamics 



(1.2.3) W(Xx, 0) = W(x, 0) , VA > 0. 

Then an hypothesis of unicity for weak solution of problem (1.2.1)-(1.2.2) shows 
that the solution W(» , •) must be self-similar. This property can be expressed 
by the relation 

(1.2.4) W(Xx, Xt) = W(x,t) , VA > 0. 

The relation (1.2.4) claims that solution W(x, t) must be searched under a 
selfsimilar form, i.e. under the form of a function of the variable j (see more 
details e.g. in the book of Landau and Lifchitz [LL54]) : 

(1.2.5) W(x,t) = U(£) , £ = |. 

• As a first step, we consider regular solutions £ i — > U(£) of the Rie- 
mann problem (1.2.1) (1.2.2). We introduce the representation (1.2.5) inside 
conservation law (1.2.1) and we obtain by this way : 

(1.2.6) dF(U(£)).^ = 

We deduce from (1.2.6) that we are necessarily in one of the two following cases : 
either vector ^ is equal to zero or this vector is not equal to zero. In the 
first case, the solution is a constant state and in the second opportunity, the 
vector ^ is necessarily equal to some eigenvector R(W) of the jacobian matrix 

dF(W). In this second case, we have the classical relation between the jacobian 
matrix, eigenvector and eigenvalue X(W) : 

(1.2.7) ' dF(W).R(W) = X(W)R{W) 

in the particular case where W = U(£). By identification between the two 
relations (1.2.6) and (1.2.7), we deduce that the vectors and R(U) are 
proportional and we deduce also : 

(1.2.8) X(U(0) = £. 

We have derived the conditions that characterize a so-called rarefaction wave. 

• It is also possible to suppose that function £ i — > U(£) admits some 
point of discontinuity at £ = a. It satisfies necessarily the Rankine-Hugoniot 
relations that link the jump of state [U], the jump of the flux [F(U) ] and the 
celerity a of the discontinuity profile : 

(1.2.9) [F(U)] = a[U]. 

Recall that the Rankine-Hugoniot relations express that the discontinuous func- 
tion £ i — > U(£) is a weak solution of conservation (1.2.1) (see e.g. the book of 
Godlewski and Raviart [GR96]). By this way, we are deriving a so-called shock 
wave. 

• We have observed in the previous section that the construction of a rare- 
faction wave is associated with some eigenvalue of the jacobian matrix dF(W). 
In the case of the Euler equations of gas dynamics, these eigenvalues are simply 
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computed with the help of the so-called nonconservative form of the equations 
that are obtained with the introduction of specific entropy s (see [LL54] for the 
proof) : 



(1.2.10) 
(1.2.11) 

(1.2.12) 
Then we set 
(1.2.13) 



dp 

dt 

du 

~di 
ds 

dt 



+ u 



dp 

dx 
du 



+ P 



du 



— + u 7T + -7^ = 
ox 

ds 

777 + U — 



dx 
1 dp 

p dx 







dx 



Z{W) = (p,u,s)' 



I u 



(1.2.14) B(Z) 



= 







c 



\ 



u 



s) 



1 dp 
- -ttKP, 

p OS 

u 

and the Euler equations that took the form (1.2.10)-(1.2.12) can be also written 
with the above matrix : 



/ 



(1.2.15) 



dz 



+ B(Z) 



dz 



= o 



dt dx 
Celerity c for sound waves used in formula (1.2.14) is defined by 



dp 



(1.2.16) 

and for a perfect polytropic gas with ratio 7 of specific heats, we have 
(1.2.17) 




c = 



• If we wish to diagonalize the matrix dF(W), it is sufficient to make this 

job for matrix B(Z) because they are conjugate as we show in the following. 
We introduce the vector variable Z inside equation (1.2.1) : 

dZ dZ 

(1.2.18) dW(Z).— + dF(W).dW(Z).— = 

and by comparison with relation (1.2.15) we have necessarily : 

(1.2.19) B(Z) = {dW{Z))- 1 .dF(W).dW{Z). 

This conjugation relation shows that R(Z) = {dW {Z))~ x • R{W {Z)) is an eigen- 
vector for the matrix B(Z(W)) with the eigenvalue X(W) if R(W) is some 
eigenvector of the jacobian matrix dF(W) that satisfies the relation (1.2.7). 
The diagonalization of matrix B(Z) is straightforward. We find 

(1.2.20) Ai = u - c < A 2 = u < A 3 = w + c 
with associated eigenvectors given by the following formulae : 
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(1.2.21) 




We remark that the derivation of scalar field Ai (respectively A3) inside direction 
Ri (respectively R3) is never null 

(1.2.22) dAi(W(Z)) • Ri(W) ^ 0, d\ 3 (W(Z)) . R 3 (W) ^ 0, VW 
whereas we have the opposite situation for the second eigenvalue A2 : 

(1.2.23) d\ 2 (W(Z)).R 2 (W) = 0, VW . 

For this reason, we will say that the eigenvalues Ai and A3 define genuinely 
nonlinear fields whereas the eigenvalue A2 = u defines a linearly degen- 
erate field. Rarefaction and shock waves are always associated with genuinely 
nonlinear fields and in what follows, we distinguish between 1-rarefaction wave 
and 3-rarefaction wave. 

© 

r i( w o)a J!, 

U 1^0>= W 

Figure 1.1 Rarefaction wave associated with eigenvalue \\ = u — c in the 
space of states ; the curve is everywhere tangent to eigenvector R\(W). 
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*l(W 1 ) = § 1 



X 



Figure 1.2 Rarefaction wave U\{£) associated with eigenvalue \\ = u — c 
inside the space-time plane. 

• A 1-rarefaction wave is a function £ \ — > Ui(£) that satisfies 

(1.2.24) — (C7i(^)) proportional to ^1(^1(0) 

and we have an analogous definition for a 3-wave £ 1 — ^ ^3(0 : a? (^3(0) * s 

proportional to -R3 (Us(£)) . Then if we integrate the vector field R\{Z) i.e. if 

we solve the ordinary differential equation ^ = Ri(Z{£)) , we find a function 

£ 1 — > Z(Ui(^)) that defines a 1-rarefaction wave, say for variable £ contained 
between two limiting values £0 and £1. We can draw the solution of differential 
equation (1.2.24) inside the space of states W ; we find a curve £ 1 — > Ui(£) 
that satisfies equation (1.2.24) and the initial condition 

(1.2.25) Uifa) = W . 

Moreover this curve is defined up to £ = £1 where the final state W\ is achieved 
(see Figure 1.1) 

(1.2.26) U^Ci) = Wl 

In space-time (x, t), we observe that we have to consider three rates of flow for 
celerities £ = ~ : 

(1.2.27) Ui(C) = W for £ < £0 (constant state) 

(1.2.28) Ui(£) variable for £0 < £ < £1 (rarefaction wave) 

(1.2.29) UiH) = W 1 for £ > £1 (constant state) 

as proposed at Figure 1.2. Notice that relation (1.2.8) imposes a particular value 
£0 for celerity : 

(1.2.30) £ = u(Wo)-c(W ) 

and an analogous value for end-point £1. It is simple to show (see e.g. Smoller 
[Sm83] or Godlewski-Raviart [GR96]) that inequality 
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(1.2.31) £ < 6 

is necessary if W\ is linked to Wq through a -1- or a 3-rarefaction wave. In 
the particular case of a 3-rarefaction wave, relations (1.2.27) to (1.2.29) are still 
correct, but relation (1.2.30) must be replaced by 

(1.2.32) £ = u(W ) + c(W ) 

and we have also an analogous relation for £i . 



1.3 Riemann invariants and rarefaction waves. 

• The practical computation of the curve £ i — > Ui(£) satisfying the 

relations (1.2.24) and (1.2.25) uses the notion of Riemann invariant. By defini- 
tion, a 1-Riemann invariant (respectively a 3- Riemann invariant) is a function 
W i — > /3 1 (H / ) (respectively W i — > (5 s (W)) that is constant along the curves 
of 1-rarefactions (respectively along the curves of 3-rarefactions) and satisfies 
by definition : 

(1.3.1) dfi(W).R 1 (W) = VW 

(respectively d/3 3 {W) . R 3 (W) = for each state W). If we express this relation 
in terms of nonconservative variables Z, we set /3(Z) = f3{W(Z)) and we have : 

d/3(Z).R(Z) = d/3(W(Z)).dW(Z). (dW(Z)y 1 . R(W(Z)) = if f3(.) is 
a Riemann invariant for the field associated with eigenvector R(»). We deduce 
from the previous calculus the relation : 

(1.3.2) d/3 1 (W(Z)).^ 1 (Z) = VZ. 

Taking into account the particular form (1.2.21) of the vector Ri(Z), the two 
following functions 

(1.3.3) p\(W) = s 

(1.3.4) fil(W) = u+ f P C -^d9 = u + — — r 

Jpo ^ 7-1 

are particular 1-Riemann invariants ; they satisfy together the relation (1.3.2). 
The expression of the 3-Riemann invariants is obtained by an analogous way : 

(1.3.5) P\{W) = s 

(1.3.6) Pl(W) = u- f C -^d0 = u - 

Jpo 7-1 

The states W on a 1-rarefaction wave issued from state Wq are explicited by 
using the relation (1.3.1) for the Riemann invariants proposed at relations (1.3.3) 
and (1.3.4). This fact express that on a 1-rarefaction wave, the two associated 
Riemann invariants are constant ; we have 

(1.3.7) s = s 
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(1.3.8) 
(1.3.9) 



u + 

z = 



2c 



U + 



2 c 



7— 1 " ' 7— 1 
u — c , £ > wo — Co . 

• We detail now the particular algebraic form that takes the description of 
the link beween a state W and its initial datum Wq through a 1-rarefaction 
wave inside a perfect polytropic gas. We first set (we refer for the details to the 
book of Courant and Friedrichs [CF48]) that inside a 1-rarefaction wave, the 
pressure p is a decreasing function of velocity : 

(1.3.10) p < po , u > uq , W issued from Wq mo a 1-rarefaction 
and due to the expression of the entropy s for a polytropic perfect gas, 

(1.3.11) s is a function of variable — 

it comes, taking into account the relations (1.2.17), (1.3.6) and (1.3.7), 

\_ 

7-1 



(1.3.12) u - u + 



(1-rarefaction) 



P Al ~ Vq 

parameterized by the non-dimensional coefficient /x > [note that this param- 
eter fi has nothing to do with the chemical potential, even if the same letter is 
used !] defined by 

2 7-1 



(1.3.13) ^ 



7 + 1 



In the plane (w, p) composed by velocity and pressure, the graphic representation 
of the 1-rarefaction (relation (1.3.12) under condition (1.3.10)) is proposed on 
Figure 1.3. 
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• For a 3-rarefaction wave, we look for (due to reasons that will be ex- 

plained in what follows) an upstream state W linked to a downstream state 
W . We obtain 

(1.3.14) s = s 

(1.3.15) u — = Uq — 

7 — 1 7 — 1 

(1.3.16) £ = w + c, £<wo + co- 

It is also easy to prove that in a 3-rarefaction, velocity is a nondecreasing 
function of pressure, i.e. : 

(1.3.17) p < po , u < Uq , Wo issued from W via a 3-rarefaction. 

A computation analogous to the one presented for 1-rarefactions shows 

j_ 

-u 4 p 21 ( ¥± 

u — uo — - — 5 — I p 2 ~f — p 7 J = (3-rarefaction). 

The comparison between relations (1.3.12) and (1.3.18) induces us to set 



(1.3.18) 



1 

— I ^7 
Po 



' 7-1 



Po 



7-1 
27 



(1.3.19) v(p; Po , Po ; 7) ; 

The graph in the plane (u, p) of the curve of 3-rarefaction (equation (1.3.18) 
under the constraint (1.3.17)) is presented on Figure 1.4. 






Figure 1.4 Curve in the plane of velocity and pressure showing the set of 
states W such that the 3-rarefaction wave that begins at state W ends at the 
particular state Wo- 
rn It can be usefull to precise the (varying) state inside a 1-rarefaction 
wave as a function of celerity £ = j . Taking into account the relation (1.2.8), 
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we have on one hand : 
(1.3.20) u - c = f 

and on the other hand taking into account the relation (1.3.8), it comes easily, 

for u\ — ci < £ < U\ — C\ : 

7-1 2 2 

1.3.21 u = ±— Ul + —ci + — £ 

7 + 1 7+1 7 + 1 

-7-1 2 7—1 

(1.3.22) c = I— m + — c, - 

7 + 1 7+1 7+1 

The elimination of pressure among the relations (1.2.17) and (1.3.7) allows the 
evaluation of density. It is the same set of operations for a 3-rarefaction wave. 
The relation (1.2.8) gives the expression of the celerity of the 3- wave as a function 
of the data : 

(1.3.23) u + c = £ 

and the Riemann invariant (1.3.15) allows the evaluation of velocity u and sound 
celerity c as a function of the data : 

(1.3.24) u = ^u r - -^-c r + -4t£ 

7 + 1 7+1 7+1 

7-1 2 7-1 

(1.3.25) C = -^—Ur + — -Cr + ^—T £ 

7 + 1 7+1 7+1 

and under the conditions U2 + C2 < £ < u r + c r . 

1.4 Rankine-Hugoniot jump relations and shock waves. 

• We consider now two particular states Wq and W linked together via the 
Rankine-Hugoniot jump relations (1.2.9) and the associated entropy inequality. 
We first remark that the Galilean invariance of the equations of gas dynamics 
allows us to think the Physics inside the reference frame with a velocity exactly 
equal to the celerity a of the discontinuity. Then the jump equations (1.2.9) 
are written again under the more detailed form 
(1.4.1) 
(1.4.2) 

(1.4.3) p(u - a) \ i ■ "' y ) + p(u - a) = 



p(u-a)} = 
p {u - a) 2 + p] = 

p{u-a)ye + J + p (u - a) 



as it is also easy to derive directly. We consider the classical expression for mass 
flux that cross the shock wave : 

(1.4.4) m = p(u- a) . 

The associated conditions of increasing physical entropy through a shock wave 
(see e.g. [CF48] or [LL54]) claim that we have 

(1.4.5) m > through a 1-shock 

(1.4.6) m < through a 3-shock. 
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The numbering of shock waves can be explained by an argument of continuity 
as follows. If the jumps inside relations (1.4.1) to (1.4.3) are weak enough, it is 
possible to show (see details in [Sm83] or [GR96] for example) that the celerity 
a of the discontinuity relative to both states Wo and W has a limit value 
equal to the common value u — c for a 1-shock and u + c for a 3-shock wave. 
The particular case where the mass flux m is equal to zero will be considered 
afterwards and corresponds to a slip surface or contact discontinuity. 

• We detail the algebra that is necessary in order to express that the state 
W is obtained from an upstream state Wq through a 1-shock wave with celerity 
a. We first remark that the entropy condition implies the following family of 
inequalities : 

(1.4.7) p > p , p > p , u-c<a<u -c , s > s 
when W is issued from Wo through a 1-shock wave. We denote also by h and 
r the specific enthalpy and the specific volume 

(1.4.8) h = e + - 

P 

(1.4.9) r = -. 

P 

From relations (1.4.1) to (1.4.3) we deduce, taking into account (1.4.4), (1.4.8) 
and (1.4.9) : 

(1.4.10) [u] = m [t] 

(1.4.11) m 2 [t] + [p] = 

(1.4.12) [h] +^[r 2 ] = 0. 

In the particular case of a polytropic perfect gas, we eliminate the mass flow m 
from relations (1.4.11) and (1.4.12) and we express the enthalpy as a function 
of pressure and specific volume, i.e. 

(1.4.13) h = -^-rpr. 

7-1 

After some lines of elementary algebra we express the specific volume r down- 
stream to the shock as a function of upstream data, downstream pressure and 
variable p introduced at relation (1.3.13). We obtain : 

(1.4.14) r = — ■ 5 — T . 

p + p z p 

We report this particular expression inside relation (1.4.11) and obtain by this 
way the square of the mass flux across the shock : 

_|_ 2 

(1.4.15) m 2 = ^L^ , W issued from Wo by a 1-shock. 

(l-p 2 )T 
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Figure 1.5 Curve showing (in the velocity-pressure plane) the set of states 
W issued from the particular state Wq through a 1-shock wave. 







(3) shock / 








w • — — 










Figure 1.6 Set of states W such that the 3-shock wave links state W to the 
particular state Wq. 

• For a 1-shock wave, the relation (1.4.2) joined with (1.4.10) allows us 

to precise the jump of velocity as a function of pressure and upstream state : 

(1.4.16) u - u + J — t—^\ — r (p-Po) = 

V Po [P + P 2 Po) 
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when W is issued from Wo through a 1-shock wave. In the plane of velocity- 
pressure variables, we can propose (see Figure 1.5) the curve characterized by 
the equation (1.4.16) and the inequalities (1.4.7). 

• For a 3-shock wave, the relations (1.4.10) to (1.4.13) are still valid. We 
do not consider in what follows a state W issued from Wo via a 3-shock wave 
but we reverse the roles and consider the set of states W that are upstream 
to the particular state Wo through a 3-shock wave. We reverse the zero index 
in previous formulae. With these particular hypotheses concerning our new 
notations, a simple calculus shows that relation (1.4.15) remains still correct 
and we have 

(1.4.17) m 2 = 7 ~5T~~ , Wo issued from W by a 3-shock. 

The algebraic relation that express the jump of velocity as a function of down- 
stream state Wo and upstream pressure takes now the form 

(1.4.18) u-uo- J ( l ~J l2 2 J p-Po) = 

V po {P + P z Po) 

when Wo is issued from W through a 3-shock wave. Taking into account these 
new notations, the inequalities for entropy condition can be written as 

(1.4.19) po < p , Po < P , wo + c <cr<w + c, So < s 
if Wo is issued from W through a 3-shock wave. A 3-shock wave is graphically 
represented represented at Figure 1.6. By comparison between the relations 
(1.4.16) and (1.4.18), it is natural to define 

l~ 1-P 2 

(1.4.20) (p(p; po, Po; 7) = \ — , . 2 — r (p ~ Po) 

y po [p + p z po) 

where parameters \i and 7 are linked by relation (1.3.13). 
1.5 Contact discontinuities. 

• When we have studied the rarefaction waves, we have introduced the 
notion of linearly degenerated field at relation (1.2.23) and we have observed 
that the second characteristic field of the Euler equations of gas dynamics is 
effectively linearly degenerated. This fact means that the second eigenvalue A2 
is also a Riemann invariant : 

(1.5.1) f3 2 (W) = u. 

Taking into account the particular expression (1.2.21) of eigenvectors, it is easy 
to see that the pressure is also a 2-Riemann invariant independent from the first 
one : 

(1.5.2) P 2 (W) = p. 
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We can search a selfsimilar regular wave £ i — > U(£) solution of the differential 
equation that expresses proportionality between and R 2 (U), but in that 

case the necessary condition (1.2.8) implies that variable £ is not allowed to get 
any variation ! We have in fact by derivation of relation (1.2.8) 

(1.5.3) d\ 2 (U).^ = 1 

whereas this last quantity is identically equal to zero (relation (1.2.23)) if we 
consider the hypothesis of linearly degeneration of second characteristic field 
A 2 = u. 



s u 0= u 



W 



o 




X , t 



W 







Figure 1.7 Contact discontinuity between the two states Wq and W. 



• We remark also that if the function £ i — > U(£) is equal to some integral 
curve of vector field R 2 , i-e. 

(1.5.4) ^ = R 2 (U) 

then we have between states Wo and W the following calculus : 

f€ dU f€ 

F{W) - F(W ) = / dF(U(ri)). — drj = / dF(U(ji)) . R 2 (U(rj)) drj 

J to J to 



X 2 (U(rj))R 2 (U(rj))dr 1 



a 2 / MU(ri))&n = *2(w-w ) 



The fourth step in the previous relations is a consequence of the linear degeneres- 
cence of velocity u which is a constant along the 2-curve defined at relation 
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(1.5.9) U 2 (0 = 



(1.5.4) . We have just shown that we have a Rankine-Hugoniot jump relation 
between states Wo and W : 

(1.5.5) [F(W)] = A 2 [W] . 

In the space of states, the rarefaction waves denned by the differential relation 
(1.5.4) and the shock curves that we founded at relation (1.5.5) are identical ! 
The second characteristic field is linearly degenerate and defines a curve in the 
space of states that joins the two states Wo and W according to the relations 

(1.5.6) u = uo 

(1.5.7) p = Po . 

In space-time variables, these two states are linked through a contact disconti- 
nuity with celerity a equal to A 2 due to relation (1.5.5) : 

(1.5.8) a = u = Uq . 

Such a wave is called contact discontinuity or slip line. We have 

W , i < Uq = a 
W , £ > w = o . 
A representation in space-time of relation (1.5.9) is proposed at Figure 1.7. In 
the velocity-pressure plane, the (nonlinear !) projections of states Wo and W 
coincide, as previously established in relations (1.5.6) and (1.5.7). 

1.6 Practical solution of the Riemann problem. 

• We propose to solve the Riemann problem (1.2.1)-(1.2.2) between two 
states Wi and W r , and we remark that the general theory proposed by Lax 
(see e.g. [Lax73]) can be applied for gas dynamics. We search two intermediate 
states Wi and W 2 such that 

(1.6.1) Wi is issued from state W\ by a 1-wave 

(1.6.2) W 2 is issued from state Wi by a 2-wave 

(1.6.3) W r is issued from state W 2 by a 3- wave. 

As a first step, we restrict ourselves to the search of velocity and pressure which 
is common to both states W\ and W 2 due to the fact that the 2-wave is a 
contact discontinuity : 

(1.6.4) u\ = U2 = u* 

(1.6.5) pi = p 2 = p* ■ 

• The relation (1.6.1) means that state W\ is issued from state W\ through 
a 1-rarefaction wave (relation (1.3.10)) or a 1-shock wave (inequalities (1.4.7)). 
We have in consequence : 

ui - ui + ^{pi ; pi , pi ; 7) = , V\ < Vi 
ui - ui + (p(pi ; pi , pi ; 7) = , pi > pi . 
In a similar way, we remark that relation (1.6.3) expresses that state W 2 is the 
upstream state if a 3- wave whose downstream state is exactly the datum W r . 



(1.6.6) 
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We can use a 3-rarefaction wave between W 2 and W r (relation (1.3.17)) or a 
3- shock wave (relation (1.4.19)). It comes : 



(1.6.7) 



u 2 
u 2 



u r 
u r 



^{P2 ; Pr , Pr ; 7) = , P2 < Pr 

<P(P2 ; Pr , Pr ; 7) = , P2 > Pr 



It is sufficient to write the equation (1.6.4) that links velocities u\ and u 2 {i.e. 
ui = u 2 = u*) under the condition (1.6.5) related to the common pressures p\ 
and P2 {pi = p 2 — P*) to set completely the problem (see also Figure 1.8). 



(T) shock 



(3) shock 




(7) rarefaction 



Figure 1.8 Solution of the Riemann problem R(W\ e { t , VF r i g ht) in the velocity- 
pressure plane. 

• The numerical resolution of problem (1.6.4)-(1.6.7) can be done by New- 
ton iterations indexed by some integer k and presented on Figure 1.9. Starting 
from a given pressure pk , we easily compute with relations (1.6.6) and (1.6.7) 
velocities and u^^k respectively associated with the 1 -wave and the 3- 

wave. We evaluate the locus of intersection of the two tangent lines issued from 
the two corresponding velocities in order to define a new value pk+i for pressure 
at iteration + The initialization of the algorithm can be obtained by the 
intersection of the two rarefaction waves, i.e. by solving the following system 
with unknowns («o, Po) '■ 

(16 8) [ u - ui + ^{p ; pi, Pi; 7) = 
\ u - u r + i/){p ; p r , Pr ; 7) = . 
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Figure 1.9 Newton iterations for the resolution of the Riemann problem in 
the plane of velocity and pressure. 

Note that in relations (1.6.8), the index "O" stands for "Osher" because solving 
equations (1.6.8) is essentially what is to be done for the computation of the 
Osher [Os81] flux decomposition (see e.g. [Du87]). This system of equations, 
even if it is a nonlinear one, can be solved exactly with some explicit algebra 
and pressure po is finally evaluated thanks to the relation : 



(1.6.9) = 2 K ' —^r > 

where sound celerities q and Cy are evaluated with relation (1.2.17). We re- 
mark that relation (1.6.9) defines effectively a positive pressure if the following 
relation of non vacuum appearance is satisfied : 

(1.6.10) U r — Ul < -^—{Ci+Cr). 

7-1 
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1 -rarefaction 



vacuum 








3 -rarefaction 



u 



Figure 1.10 

plane. 



Apparition of cavitation phenomenon in the velocity-pressure 



1 -rarefaction 



3 -rarefaction 




Figure 1.11 Apparition of cavitation phenomenon. Solution of the Riemann 
problem in the space-time plane. 

• If in contrary relation (1.6.10) does not hold, the two curves associated 

with both rarefaction waves do not intersect in the velocity-pressure plane. A 
vacuum appears (see Figures 1.10 and 1.11). The solution of the Riemann pro- 
blem is no longer mathematically well defined in the sense of Lax and contains 

a zone without any matter defined by celerities £ such that 

2 ci 2c r 

1.6.11 ui + l - < £ < u r r -. 

7— 1 7— 1 

For such celerities, the pressure and the density are null whereas velocity is not 
defined. This cavitation process remains an exception but can be still completely 
solved as we have just seen. 
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• When relation (1.6.10) is satisfied, the Newton algorithm illustrated at 

Figure 1.9 is convergent towards a pair (it*, p*) composed by the common ve- 
locity it* G IR and pressure p* > of the two intermediate states W\ and 
W<i. The calculus of the density of these two intermediate states depends on the 
choice of the wave effectively used for the resolution of the Riemann problem. If 
pressure p* is less or equal to the left pressure pi, the 1-wave is a rarefaction 
wave ; then the entropy remains constant and we have 



*x 1 

P \~ 



(1.6.12) pi = ( — ) pi , f < Pi- 

On the opposite case, we use a 1-shock wave and taking into account the relation 
(1.4.9), we get finally 

n a ion P* + ^ 2 P l 

1.6.13 pi = ■ Pi ' P > Pi ■ 

pi + /i 2 p* 

The celerity <ti of the shock wave can be explicited : 



(1.6.14) cri = ui - J — — p* > pi 

" (1-P 2 ) Pz 



• For the 3-wave, there is an analogous discussion that conducts finally 

to the following relations : 



* 1 

(1.6.15) p 3 = [ — ) Pr , P* < Pr 

\p r y 



* 2 

f-i a -tc\ P ' A* Pr *^ 
(1.6.16) p 3 = ■ — Pr > P > P 



r 



/ -|- 11 Pj- 

(1.6.17) a 3 = u r + J— — - p* > p r ■ 

The qualitative comportment in space-time of the solution of a Riemann problem 
for gas dynamics is illustrated on Figure 1.12. 



23 



Francois Dubois 




2) Partial Riemann problem for hyperbolic systems. 

• In this section we generalize in two directions the notion of Riemann 
problem presented in the previous section for the particular case of gas dynamics. 
First we consider a general nonlinear system of conservation laws and second we 
introduce the notion of partial Riemann problem between a state and a 
manifold. 

2.1 Simple waves for an hyperbolic system of conservation laws. 

• We study a system of conservation laws in one space dimension. The 
unknown function IR x [0, +oo[ 3 (x, t) i — > W(x, t) G Q takes its values 
inside a convex open cone O, included in IR m (m is a fixed positive integer and 
m = 3 in the particular case of the Euler equations of gas dynamics studied 
previously) : 

(2.1.1) VWeQ, VAX), XWeQ, Q c IR m . 

The flux function 3 W i — > F(W) G IR m is supposed to be sufficiently regular 
(of C 2 class typically). The conservation law takes the following classical form 
(see e.g. Godlewski-Raviart [GR96]) : 

(2.1.2) —W(x, t) + t)) = 0. 

We suppose that the system of conservation laws (2.1.2) is a strictly hyperbolic 
system, that is for each W G £1 , there exists m real distinct eigenvalues \j(W) 
satisfying conventionnaly the ordering condition 
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(2.1.3) Ai(W0 < \ 2 {W) < ••• < X m (W), Wen 
and associated with eigenvectors Rj (W) G IR m : 

(2.1.4) dF(W).Rj(W) = Xj(W)Rj(W), Wen, j = l,---,ra. 
Each vector r G H m can be decomposed in the basis of eigenvectors Rj(W) and 
the coordinates in this basis define the family of left-eigenvectors (lj(W)) m 
which is exactly the dual basis of system (Rj(W)) ._ x m : 

k=m ^ ' 

(2.1.5) r = J2 > ') Rj(W) 

k=i 

with the classical property : 

(2.1.6) (lj(W), R k {W)) = 5 jk . 

Hypothesis 1. Genuinely nonlinear or linearly degenerate fields. 

We restrict ourselves to systems of conservation laws such that for each integer 
j with 1 < j < m, the j-th field is supposed to be 

• either genuinely nonlinear, that is 

(2.1.7) d\j(W) • Rj (W) = 1 , VW eft (j-th field genuinely nonlinear) 

• either linearly degenerate, id est 

(2.1.8) d\j(W).Rj(W) = , VW en (j-th field linearly degenerate). 

• When the j-th field is genuinely nonlinear, it is possible to construct 
the so-called j— wave issued from a particular state Wo G n. This (genuinely) 
nonlinear wave can be considered from two points of view. On the first hand, 
it is a curve e i — > % . (e ; Wo) inside the space n of all the states and on the 
other hand for real variable e fixed sufficiently small, it is possible to construct 
a self-similar weak solution of the conservation law (2.1.2) between state Wo 
and state % . (e ; Wo). 

• When e > this particular wave is a j— rarefaction wave and is defined 
as an integral curve of the vector field n 3 W i — > Rj(W) G H m in the space 
of states : 

(2.1.9) Wo)) = Rj(x.(e; W )) , e >0 

(2.1.10) X C (0,Vo) = W . 

In space-time plane, this j— curve allows to construct a continuous rarefaction, 
which is a particular solution of the conservation law (2.1.2) : 

(2.1.11) W(x,t) = Wo if y < Xj(W ) 

(2.1.12) W(x,t) = x.(^-*j(Wo);Wo) if Xj(W ) <j< Xj( X . (e ; W )) 

(2.1.13) W(x, t) = x.(e; W ) if - > ^j{x.(e; W )j . 
Moreover we have 

(2.1.14) Xj(W(x,t)) = - t if XjiWo) < - t < Xjix.ie; W )) . 
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• When e < the j-th nonlinear wave is a j— shock wave of celerity 
<7j(e; Wo) satisfying the entropy condition (e.g. in the sence of Lax [Lax73]). 
The states W and y.(e; W ) are linked with celerity Oj(e\ W ) in state space 
according to the Ranlane-Hugoniot jump conditions : 

(2.1.15) F( X .(e- Wo)) - F(W ) = a 3 (e ; W ) ( X . (e ; W ) - W ) , e<0. 

In space-time space, this discontinuous self-similar j — shock wave IR x [0, +oo[ 
3 (x, t) i — > W(x, t) E Q of strength | e | is characterized by the following two 
conditions 

(2.1.16) W(x, t) = W if - < (Tj(e; W ) 

(2.1.17) W(x, t) = y.(e; Wo) if - > ^(e; W ) . 

• When the j-th field is linearly degenerate, the construction of the j— wave 
issued from the particular state Wo E f2 is still possible. Due to the condition 
(2.1.8), the j° eigenvalue Xj(W) is constant along the integral curve defined 
in relations (2.1.9) (2.1.10) and we have also all along this curve the following 
jump relation : 

(2.1.18) F( X .(e; W )) - F(W ) = Xj(e; W )( Xj (e; W ) - W ) , Ve. 

In space-time plane, a self-similar j— contact discontinuity can be constructed 

between states Wo and y (e ; Wo) and we have : 

j 

(2.1.19) W(x,t) = Wo if X - < \j(e;W ) 

(2.1.20) W(x,t) = x .(e; Wo) if j>X 3 (e;Wo). 

All this material is summarized in the following result (see e.g. [GR96]). 
Proposition 1. 

We suppose that the hyperbolic system of conservation laws (2.1.2) satisfies 
Hypothesis 1. Then for each state Wo E Q, and for each integer j (0 < j < 
m), there exists some vicinity Qj of in H and there exists a j— wave 
X. : Qj x Q 3 (e, 

Wq) 1 — > X.( e 'i Wo) E Q, which is a regular (of C 2 class) curve. This corre- 
sponds to a j— rarefaction when the j° field is genuinely nonlinear and e > 0, 
to an admissible discontinuity (weak shock satisfying an entropy condition) of 
conservation law (2.1.2) when the j° field is genuinely nonlinear and e < or 
to a j— contact discontinuity when the j° field is linearly degenerate. Moreover, 
the j° wave y (• ; Wo) is continuously derivable at starting point Wo and we 

3 

have : 

(2.1.21) x.(e;W ) = Wo + eR j (Wo) + 0(e 2 ) 1 eE@ 3 . 

• We remark that in space-time, a j-wave defines an entropy solution of 
hyperbolic system (2.1.1). 
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2.2 Classical Riemann problem between two states. 

• We consider the solution of the particular Cauchy problem associated 

with the conservation law (2.1.1) and the particular initial datum 

(2.2.1) W( X , ) = {^ : HI 

composed by two constant states on each side of a discontinuity located at x — 0. 
In what follows, this problem is refered as the classical Riemann problem 
and is denoted by R(Wi , W r ). We have the following Theorem due to Lax (see 
e.g. [Lax 73]). 




w u =w, 



w right =w 3 



wave 



o 


1 (?) wave 
/ 


[\J wave 


i y (3) wave 




1 x^^^^W , =W 3 

vv nght vv 


~3 


X 



Figure 2.1 Solution of the Riemann problem in the state space Q and in the 
space-time plane. 
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Theorem 1. Solution of the classical Riemann problem. 

We suppose that the hyperbolic system of conservation laws (2.1.2) satisfies 
Hypothesis 1 : for j E {1, • • • , ra}, the j° field is either genuinely nonlinear or 
linearly degenerate. Then for each Wi E 0, there exists a vicinity 3^ of state 
Wi (y E V(Wi), set of all vicinities of the particular state Wi ) such that for 
each state W r lying in 3^ {W r E y), the Riemann problem R(Wi , W r ) has a 
unique entropy solution composed by at most (ra + 1) states separated by (at 
most) ra elementary waves. 

• Figure 2.1 shows what has to be done for the resolution of the classical 
Riemann problem in the particular case ra = 3 that corresponds to the Euler 
equations of gas dynamics in one space dimension. The two intermediate states 
W 1 and W 2 are linked with data W\ and W r acording to the wave relations : 

(2.2.2) W 1 = Xl (\\Wi) 

(2.2.3) W 2 = %2 ( V W 1 ) 

(2.2.4) W r = X 3 (e 3 ;W 2 ) 

and by elimination of the two intermediate states W 1 and W 2 we obtain 
(2-2.5) X 3 (e 3 ;x 2 (e 2 ;x i (e i ;^))) = W r 

which is a set of three equations with three scalar unknowns e , e and e . 

^ 1 ' 2 3 

When parameters e . are determined (in practice for gas dynamics with the 

method presented in the previous section) each wave is acting in space-time 
space as presented on Figure 2.1. 

• The proof of Theorem 1 in the general case consists in studying the 
chaining of m elementary waves, i.e. the mapping x defined by the relations 

TR m D 6i x • • • x 6 m 3 (e i , • • • , ej = e — ► X (e) E Q 

= X m (e m ; X m _ ± (e m _ 1 5 • • • 5 Xi{\ 5 Wi) • • •)) • 
Two local properties have to be derived concerning on one hand the mapping x 
itself at the origin : 

(2.2.7) x(0) = W l 

and on the other hand the tangent vector field dx considered at the same point 
has a very simple expression : 

m 

(2.2.8) dx(o).7? = J2vjRj(Wi), v= {m,---,Vm)eJR m . 

Then the local inversion theorem proves that the equation 

(2.2.9) x(e) = W r 

has a unique solution. For the details of this proof, we refer to Godlewski and 
Raviart [GR96]. □ 



(2.2.6) 
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• In space-time domain, the entropic selfsimilar solution IRx ]0, +oo[3 

(x, t) i-)- U(j ; Wi, W r ) of the Riemann problem R(Wi, W r ) is constructed as 
follows. Let 

3 3 

w m = x (e ; w m_1 ) = w r 

i s\-rn. \ rn. ' ' 



(2.2.10) 



be the intermediate states, \i- (Wi, W r ) the smallest wave celerity of the j° 
wave and n~^(Wi,W r ) the corresponding maximal wave celerity. If the j° wave 
is a rarefaction, we have from (2.1.11) and (2.1.13) : 

(2.2.11) fiJ(W h W r ) = XjiW^ 1 ), n+{W u W r ) = Xj(W j ) 

whereas in case of a j-shock wave or j-contact discontinuity, we have due to 
(2.1.16) and (2.1.17) : 

(2.2.12) fiJ(W h W r ) = fi+(W l: W r ) = ajie.-.W^ 1 ). 
Then the selfsimilar solution U(£] Wi, W r ) satisfies 

(w° = w l , i<^{w u w r ) 



(2.2.13) U{i-W h W r ) ={ 



x 

x- It 

3 \ t 



w j , 



flj(W h W r )<i<^{W U W r ) 
A*t(Wi, W r ) < £ < Hj +1 (Wu W r ) 



[W m = W r , fl+(W h W r ) 



2.3 Boundary manifold. 

• The Riemann problem is a usefull tool to prescribe weakly a boundary 

condition for an hyperbolic system of conservation laws when a right state W r 
is supposed to be given (see Section 3.3). For physically relevant conditions 
in gas dynamics, the data at the boundary are of the type "the pressure is 
known" or "total pressure and total temperature are given" , and in consequence 
do not define explicitely a single state W r . Nevertheless, a manifold of states 
is associated with these sets of incomplete boundary data, as we will develop 
in Section 3.4. This physical situation motivates the following definition of a 
boundary manifold. 
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dO>(X) • 




Figure 2.2 Tangent vectors at the particular state X of the boundary mani- 
fold .M right- 

• We restrict ourselves to a local manifold M. r around a given state Wo- 
Consider a state Wo G f2, a vicinity f2 of this state, i.e. 

(2.3.1) Wo G , fioe V(W ) , O c Q c E m 

and an invertible regular local chart $ defined on the vicinity Qq and taking 
its values in some vicinity O C H m of in H m : 

(2.3.2) tt 9 X i— ► $(X) G 6 c H m , BeV(O). 

Mapping $ is one to one, of C 1 class ($ G ^(Qq, B)) and its inverse mapping 
$ _1 is also of C 1 class ($ _1 G C X (B, fi ))- Tne boundary manifold M r i g ht = 
A^ r is here defined locally as the set of states X in the vicinity of Qo (X E Qq) 
satisfying the equations 

(2.3.3) $i(X) = ••• = $ p (X) = 0, Iefi . 

The index p is a fixed integer such that < p < m and is the co-dimension of 
the boundary manifold Ai r . Recall that for each state X (or point X G Ad r ), 
the vector in H m has its p first coordinates equal to zero : 

(2.3.4) $(M r ) C {(0, ••• ,0, y), y G IR m_p } . 
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• A system of tangent vector fields , • • • , £ m (X)) at point X E 

M. r is obtained by lifting the tangent mapping d$(X) at point X E Qo- Let 
e k be the k° vector of the canonical basis in linear space IR m . We have by 
definition : 

(2.3.5) d$(X).£ k (X) = e k , X E M r , & > p + 1 
as illustrated on Figure 2.2. We make a new hypothesis. 

Hypothesis 2. Transversality. 

We suppose that the following family E(Wo) of m vectors : 

(2.3.6) E(W ) = (Ri(W ) , • • • , R P (W ) , £ p +i(W ) , • • • , ^ TO (W )) 
is a basis of linear space H m . 



2.4 Partial Riemann problem between a state and a manifold. 

• The partial Riemann problem P(Wi, M. r ) between a state W\ and a 

manifold A4 r is by definition the Cauchy problem for the hyperbolic system of 
conservation laws (2.1.2) associated with the following constraints relative to the 
initial condition : 

(2.4.1) w { ,, 0) HI 

The above definition of the partial Riemann problem P(Wi, A4 r ) has been first 
proposed in a particular case in [Du87] . We have also used it in [Du88] , [DLF89] , 
[DLL91] and [CDV92]. We first remark that if p = m and M r = {W r } , then 
the partial Riemann problem P{Wi, M. r ) reduces to the classical Riemann 
problem R(Wi, W r ). The following result has been first presented in [Du98]. 

Theorem 2. Existence of a solution for the partial Riemann problem. 

We suppose that the hyperbolic system of conservation laws (2.1.2) satisfies 
Hypothesis 1 and that the boundary manifold M. r is defined as above with a 
given state Wq E Q that satisfies Hypothesis 2, a vicinity Qo of Wo (^o £ 
V(W )) and a local chart ft 9 X i — > <$>(X) E 6 C H m , 6 E V(0) and 
satisfying the relations (2.3.3)- (2.3.4). Then there exists a vicinity Qi C flo 
of state Wo such that for every state Wi E fii, the partial Riemann problem 
P(Wi, M. r ) defined by (2.1.2) and (2.4.1) has an entropy solution in space Q 
composed by at most (p + 1) states 

W l = W",W 1 = Xl {\ , W°) , • • • , W p -' = x v _ x (e^ , W?- 2 ), 

W p = Xp (e p , Wv- 1 ), W?eM r , 
separated by (at most) p simple waves. In space-time, the solution Hx]0, +oo[ 
3 (x, t) i — > U(j ; Wi, M. r ) is obtained by superposition of the p simple waves 
considered previously. With the notations proposed in (2.2.11) and (2.2.12), we 
have 



(2.4.2) 
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r w° = wi , £ < 



(2.4.3) [/(£; Wi, W r ) = { 



x 

X. It 

3 V t 



r 5 



+ 




w°=w 



p=2, m=3 





(l^ wave 


(7) wave W 1 


/ w 2 

/ 

' W 2 G M • u, 
yv ^ 11 right 

f ► 


w°=w ^\ 

left X. 




X 



Figure 2.3 Resolution of the partial Riemann problem P(W\ e ft : Alright) in 
the state space £1 and in the space-time plane for a manifold M. r i s ht of codi- 
mension p = 2. 
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Figure 2.4 Resolution of the partial Riemann problem P(Wi e ft, Alright) in 
the state space Q and in the space-time plane for a manifold Alright of codi- 
mension p — 1. 

• The proof of Theorem 2 follows the construction illustrated on Figures 

2.3 and 2.4 for two particular cases. The idea is to apply the implicit function 
theorem to the mapping defined as follows. We first introduce a vicinity A = 
6i x • • • x Q p of in W : 

(2.4.4) JR p D 6i x ••• x9 p = A 3 ^ , ••• , e p ) = e, A e V(0) , 

and a sub- vicinity fl± C flo of given state Wo in such a way that the mapping 
ty p defined by the chaining of the first p simple waves is well defined : 
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f A x n, 3 (e, W) ► tf (e, W) e fi 
( 2 - 4 -5) ( *,(«, W) = Xp (e, ; x,., \r: V 'I':-")). 

• We iterate this mapping with the p first components $1, • • • , $ p of the 
local chart $ defined at relations (2.3.2) to (2.3.4) : 

J A x Oi 3 (e, W) ,— ► <p(e : W) G H p 
^• 4 "^ \^(e, W) = {^(%(e,W)),-- - ,%(%(e,W))). 
We construct a solution of the partial Riemann problem composed by the con- 
servation law (2.1.2) and the constraints (2.4.1) for the initial conditions with 
the p first waves issued from the left state W\ . In other terms, we search a right 
state W r G M. r as in the Riemann problem R(Wi, W r ) id est under the form 

(2.4.7) W r = %(e,Wi), W r eM r . 

The determination of the state W r satisfying the conditions (2.4.7) is equivalent 
to the research of parameter e G A c H p that satisfy the equation 

(2.4.8) ( p(e 1 Wi) = 0. 

• We have the following natural property : 

(2.4.9) tf p (0, W ) = W G M r , v(0, Wo) = G H p 

and moreover from Hypothesis 2, the m x m matrix composed by the family 
E(W ) = (Ri(W ) , • • • , R P (W ) , fp+i(W ) , • • • , U(Wo)) has a rank equal 
to m. Then the same property holds after applying the linear one to one mapping 
d<I>(W„) : 

(2.4.10) rank(d$(W ). J Ri(W ), ••• , d$(W ) . R P (W ) , e p+1 , • • • , e m ) = 
m 

due to the definition (2.3.5) of tangent vectors £ p+ i(Wo) , • • • , £m(Wo). In con- 
sequence, when we look to the matrix defined at the relation (2.4.10), we observe 
that the block composed by the p first lines and the p first columns at the top 
and the left of this matrix has a rank exactly equal to p. Moreover, this p x p 
matrix is exactly equal to (0, Wo) and in consequence this jacobian matrix 
is invertible : 

(2.4.11) — — (0, Wo) is an invertible pxp matrix. 
oe 

• Now the theorem of implicit functions proves, with the eventual con- 
straint that vicinity fii may have to be reduced, that the equation 

(2.4.12) <p(e, W) = 

admits a unique solution (e, W) in the vicinity of (0, Wo) and it takes the 
form e = ir(W) : 

, v ( 3 n ± 3W ^ e = tt(W) G A c M p , 7r G C(Oi, A) such that 
^•4.13) | w ^ = o, WGOx) (e = tt(W)). 

The general structure of the solution of equation (2.4.8) (or of the equivalent 
equation (2.4.12)) is a consequence of the implicit function theorem presented in 
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(2.4.13) : it allows the determination of the strengh e of the waves as a function 
of left state Wi. Then Theorem 2 is established. □ 

• The extension of the previous notion to a partial Riemann problem 

P(Mi, W r ) composed by a boundary manifold A^i e ft = Mi and a state 
Wright = W r is defined by the conservation law (2.1.2) and the initial conditions 

(2.4.3) W(x,0) [ =v l r \ x>0 . 

If codim(A^;) = p, the construction proposed at Theorem 2 can be extended 
without difficulty. For each state W r sufficiently close to the manifold Mi, the 
partial Riemann problem P(Mi, W r ) admits an entropy solution composed by 
the last p waves of the Riemann problem : 

f W°eMi, W 1 = X ,Ae , 1} W°), ■■■ , 

(2 4 1 5) < m-p+l m-p+l 

{ '' } \ W?- 1 = X Ae W p - 2 ), W? = X^^W*- 1 ) = W r . 

\ m — l m — 1 m m 



2.5 Partial Riemann problem with an half- space. 

• The notion of partial Riemann problem can be extended to the partic- 
ular situation of a "half space" defined as follows. We first introduce the set 
M.right = M r defined by 

(2.5.1) Mr = {WeQ, Xi(W) > 0}. 

The partial Riemann problem P(Wi , M r ) between the particular state Wi and 
the half space M r is still defined by the partial differential equation (2.1.2) and 
the initial constraints (2.4.1). We have the 

Proposition 2. 

Partial Riemann problem with a particular half space. 

We suppose that the 1-field is genuinely nonlinear. A particular solution of 
partial Riemann problem (2.2.1)-(2.4.1) between left state Wi and right half 
space M r defined in (2.5.1) can be constructed as follows : 
. if Ai(Wj) > 0, then W t E M r and Hx[0, +oo[^ (x, t) \ — > W(x, t) E Q 
is equal to a constant state : 

(2.5.2) W(x, t) = Wi, Ai(Wz) > 0. 

• if \\{Wi) < 0, the following rarefaction wave : 



(2.5.3) W(x, t) 



Wi if | < Ai(Wi) 



%i (| - x 1 (w l ) , wi) if \AWi) < | < o 
^ x^-x^Wi), w t ) if - t > 

is an entropy solution of the partial Riemann problem (2.1.2) (2.4.1) (2.5.1). 
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• The proof of Proposition 2 is a direct consequence of the notion of rare- 
faction wave explicited at relations (2.1.11) to (2.1.13). When A4 r is the half 
space of "supersonic outflow" (Ai = u — c > for the Euler equations of gas 
dynamics), the solution of the partial Riemann problem is trivial when state 
Wi satisfies this condition. On the contrary, the solution of the partial Riemann 
problem is constructed with the help of a 1-wave for linking the "subsonic left 
state" Wi (ui < ci) to the half space M. r . The "end-point" of the rarefaction is 
the sonic state W r = x i ( — ^i(Wi) , Wi) E M. r for particular celerity j = 



(see Figure 2.5). 



© 



w°=w 



half space 
(p=0, m=3) 



□ 












(^T^ rarefaction \ 






Wright e M right 
► 
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Figure 2.5 Resolution of the partial Riemann problem P(Wi e ft, M-right) in 
the state space £1 and in the space-time plane for an half space -M r i g ht describing 
a supersonic outflow. 
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3) Nonlinear boundary conditions for gas dynamics. 

• In this section we use the notion of partial Riemann problem to consider 
boundary conditions for gas dynamics. We first make the link between this 
notion and the way that linear hyperbolic systems are well posed in the sense 
of least squares. Then we detail some cases where physically relevant boundary 
conditions can be interpreted with particular partial Riemann problems. 

3.1 System of linearized Euler equations. 

• Most of the known mathematical results concern hyperbolic linear equa- 
tions. We linearize the Euler equations of gas dynamics around a constant state 
Wo with density po, velocity wq, pressure po and sound celerity cq. We set : 

(3.1.1) W = Wo + W 

and we neglect second order terms relatively to the variable W' . In particular, 
the incremental variables p' = p — po, u' = u — uo, s' = s — so define the 
incremental vector Z' of nonconservative variables : 

(3.1.2) Z' = (p'.u'.s'f 

and the difference of pressure p' = p — p is given at the first order as a function 
of the incremental thermodynamic variables p' and s' : 

(3.1.3) p' = p -po = clp' + j- s {Wo)s'. 

Starting from expression (1.2.15) of the Euler equations of gas dynamics, we get 
at the same level of approximation : 

(3.X.4) °* +B(Wo) °* = o. 

• As in section 1.2, we diagonalize matrix B{Wq) whose expression has 
been given at relation (1.2.14) and eigenvectors RjiWo) in (1.2.21). We can 
express the components ipj of incremental vector Z' in the basis of vectors 

Rj(Wo). These variables ipj are called the characteristic variables : 

3 

(3.1.5) Z' = ^tpjRjiWo) = tp.R(W ) 

3 = 1 

and we have from (1.2.21) the following expressions : 

Pi = 7, 2 (P' ~ POCQU') 



(3.1.6) 



2p c z Q 
1 , 

(f2 = ~~2 S 



^3 = W—ji(p' + Pocou'). 

IPO Cq 



• The change of variables IR 3 3 Z' \ — > cp G IR 3 allows to decouple the 

system of equations (3.1.4) into three uncoupled advection equations : 
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where A (Wo) = diag (uo — cq , uq , uq + cq). The system (3.1.7) is called the 
characteristic form of the linearized Euler equations. The above study motivates 
the mathematical study of linearized hyperbolic linear systems. 

3.2 Boundary problem for linear hyperbolic systems. 

• In a classical article, Kreiss [Kr70] developed the notion of well posed 
problem for initial boundary value problems in the quarter of space x < L 
and t > associated with linear hyperbolic systems of the type (3.1.4) or (3.1.7). 
We suppose that the boundary x = L is non-characteristic, id est 

(3.2.1) Wo -c ^0, w ^ 0, w + c ^0, 

and we denote by Aq~ (respectively A^) the negative part (respectively the 
positive part) of the nonsingular matrix A(Wo). We decompose also the charac- 
teristic vector (p into the ingoing components (p~ and the outgoing components 

= (p~ + (f + 

(3.2.2) {<p~ = { &;)i<j<3 , <pj = <Pj if Xj(W ) < 0, ifj = if Xj(W ) > } 
= { (y+)i<i<3 , ^ = <Pj if A,(W ) > 0, <p+ = if Xj(W ) < 0} . 




t 



(p" (m-p componants) 

V 

\ 

\ 

^^^^^^^^^ 


i 

(p + (p componants) 




L x 



Figure 3.1 Characteristic directions at the boundary x = L of the domain 
[0, L\. 

Due to our choice to consider the exterior of the domain "at the right" of the 
domain {x < L] of study, the vector tp~ is associated with negative eigenvalues 
of matrix A(Wo) and vector cp + corresponds to positive eigenvalues of the same 
matrix. With the above notations, Kreiss has proved [Kr70] that in very general 
situations, the linear hyperbolic system associated with the initial condition 
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(3.2.3) ip(x, 0) = cpo(x), x<L 

and the boundary condition parameterized by the reflection operator E that 
associates to the outgoing characteristics variables ip + a linear function E • cp + 
which is an incoming caracteristic variable : 

(3.2.4) c^- = E. x = L 
is well posed in space I? . 

• The boundary condition (3.2.4) can be interpreted in terms of charac- 
teristic directions : the field <p~ along the boundary is an affine function of 

the outgoing characteristics cp + (see Figure 3.1). The previous result remains 
true in more complicated situations. In particular, the multidimensional case 
can be considered with the same arguments except that direction x has to be 
replaced by the normal direction (see e.g. Higdon [Hi86]). In the characteris- 
tic case where e.g. the reference velocity uq is null, Majda and Osher [M075] 
have extended Kreiss' result. We remark also that the boundary condition can 
naturally be written in terms of a boundary manifold, even if it is an affine 
manifold in the present case. We introduce the incremental input vector Z~ 
and the associated output vector Z + : 

z' — z~ + Z^ 

Z~ = £ vjRjIWo) = ^-i^o) 

(3.2.5) <J Xj(w )<o j 

z+ = J2 VjRj( w o) = 2>/W°)' 

A i (W )>0 j 

and define the manifold M. r i g ht = M. r by the condition 

(3.2.6) M. r = { Z' given by relations (3.2.5), ip~ = E«(£> + + g } 
and re-write the boundary condition (3.2.4) under the equivalent form : 

(3.2.7) Z' £ M r . 

• When we consider the boundary condition (3.2.7), we distinguish classi- 
cally between four cases, following the sign of velocity uq and the modulus | uq \ 
compared with sound celerity cq. If uq < at x = L, the fluid enters inside 
the domain ] — oo, L] and if uo > the boundary {x = L} is an output. If 
| wo | < Co, the flow is subsonic whereas if | uq \ > Co, it is supersonic. 

(i) Supersonic inflow (w < — Co). The outgoing component of </?, 
i.e. ip + is reduced to zero, the manifold M. r is reduced to the unique point 
g»R(Wo) and M. r is of codimension p = 3. The boundary condition (3.2.7) 
is equivalent to prescribe all the components of vector Z' . 

(ii) Subsonic inflow (— Co < Uq < 0). The number of incoming char- 
acteristics is two and there is one outgoing characteristic direction. The manifold 
A4 r is of co-dimension p = 2. The linearized problem is well posed when any 
of the following pairs of two variables are given (see Oliger-Sundstrom [OS78] 
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or Yee-Beam- Warming [YBW82]) : (density, pressure), (velocity, pressure) or 
(enthalpy, entropy) ; we extend this context in the next section. 

(iii) Subsonic outflow (0 < Uq < Co). Only one characteristic is going 
inside the domain of study and two are going outside. The manifold M. r is of 
co-dimension 1 and only one relation has to be given between the scalar data 
on the boundary ; we remark that the classical choice of imposing the pressure 
p = Po can be written after linearization : 

(3.2.8) p' = 0. 

This condition is equivalent to choose 

(3.2.9) £ = (0, -1), g = 

inside relations (3.2.5) and (3.2.6) because the boundary condition (3.2.8) can 
also be written under the equivalent form : 

(3.2.10) if! = -tp 3 . 

(iv) Supersonic outflow (uo > Cq). All the characteristic directions 
are going in the direction opposite to the domain {x < L} (id est cp~ = 0) and 
condition (3.2.7) does not carry any information : no numerical datum has 
to be prescribed for a supersonic outflow. 

• As a complement of the previous cases, we can add the important case 
where uq = that corresponds physically to a rigid boundary. The linearized 
rigid wall boundary condition takes the form : 

(3.2.11) u' = 

and the latter boundary condition is equivalent to prescribe a null mass flux 
through the boundary for the linearized equations. Note that the previous con- 
dition (3.2.11) can also be written on the form (3.2.4) with the particular choice : 

(3.2.12) E = (0, 1), g = 

and conditions (3.2.11) (3.2.12) correspond to the relation 

(3.2.13) i Pl = <p 3 

between the characteristic variables. It can be shown (e.g. Oliger-Sundstrom 
[OS 78]) that the the initial boundary value problem (3.1.4), (3.2.3) and (3.2.12) 
is well posed with the natural L 2 condition. 

3.3 Weak Dirichlet nonlinear boundary condition. 

• For two particular systems of conservation laws of the type 

dW d 

^ -ar + 2*^ = 0, 

i.e. for linear hyperbolic systems and nonlinear scalar conservation laws, fol- 
lowing ideas developed by T. P. Liu ([Li77], [Li82]) for initial boundary value 
problems and transonic flow in nozzles, we have remarked with Le Floch in 
[DLF87] and [DLF88] that an efficient way to set a boundary condition of the 
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type "given state W r i g ht outside the domain of study" , say {x < L} to fix the 
ideas, is to consider a weak form of the Dirichlet boundary condition W = W r 
in the following sense : 
(3.3.2) W(L~,t) E B{W r ). 

In relation (3.3.2), the state W(L~, t) is the limit value of internal state W(x, t) 
as x < L tends to the boundary, i.e. x = L~ and B(W r ) is a set of ad- 
missible states at the boundary associated with datum Wo- Even if other 
formulations are possible and are still in development (see among others Nishida- 
Smoller [NS77], Bardos, Leroux and Nedelec [BLN79], Audounet [Au84], Benab- 
dallah [Be86], Benabdallah-Serre [BS87], [DLF88], Bourdel, Delorme and Mazet 
[BDM89], Gisclon [Gi94], Serre [Se96]), we focus here on the use of the Riemann 
problem for taking into account a physical boundary conditon, in particuler for 
gas dynamics. We define the set of admissible states at the boundary as follows : 

q q\ r(w \ / vauies a t f = 0~ °f the entropic solution 1 
[6.6.6) &{w r ) - | of the Riemann problem R ^ Wl Wr ^ wen]' 

• The advantage of the previous definition is that the initial boundary value 
problem parameterized by constant boundary datum W r , constant initial 
condition Wo and set on space-time domain (x, t) E] — oo, L[x[0, +oo[ by the 
conditions 

dW d 

+ — F(W) =0 x < L, t>0 



(3.3.4) 



Ot Ox 

W(x, 0) = W x < L , t = 

W(L~, t) e B(W r ) t > 

is well posed in conditions analogous to Theorem 1 for the Riemann problem 
in Lax's theory (see Section 2.2). Moreover, the set B(W r ) is easy to evaluate 
exphcitely. In [DLF88], we have calculated the boundary set B(W r ) in the 
particular case of Euler-Saint Venant equations of isentropic gas dynamics. Even 
when the state W r corresponds to a supersonic inflow in the linearized analysis, 
the admissible set B(W r ) is reduced to {W r } only in a vicinity of state W r . For 
large differences between limit state W(L~, t) and weakly imposed boundary 
state W r these two states can be linked together with an entire family of waves 
that compose the Riemann problem, and state W(L~, t) can even correspond 
to a supersonic outflow ! 

• We propose here to extend the previous weak boundary Dirichlet con- 

dition (3.3.4) to a manifold M. r . We first introduce a new set (3(M. r ) of 
admissible states at the boundary : 

(3 3 5) B(A4 ) { vames a ^ f = 0~ °f an entropic solution of the 1 
^ ' * ' ^ r > ~ \ partial Riemann problem P(W, M r ), W E CI J 

We remark that this set of admissible states is a natural extension of definition 
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(3.3.3) relative to admissible states associated with a single state W\. The ini- 
tial boundary value problem associated with datum M. r is now formulated as 
follows : 



(3.3.6) 



OW d 

-dt + Tx F{w) = X<L > t>0 

W(x, 0) = W x < L , t = 

W(L~, t) £ P{M r ) t > 0. 



• Theorem 2 shows that when constant initial datum Wo and boundary 

manifold M. r are closed enough, the problem (3.3.6) is well posed in the family 
of solutions composed by the p "first" waves as described in relation (2.4.2). The 
solution of the initial boundary value problem (3.3.6) is the self-similar solution 



(x, t) i — y U{^ k ; Wi, Mr) described at relations (2.4.3) in Theorem 2. We 
remark that only the waves with negative celerities contribute to problem (3.3.6) 
even if the partial Riemann problem contains waves with positive celerities (see 
e.g. figure 2.3 or 2.4). 



3.4 Some fluid boundary conditions. 

• We show in this section that for a family of classical examples, it is possi- 
ble to introduce a partial Riemann problem in order to consider fluid boundary 
conditions for the gas dynamics equations. We examine five particular cases : 
given state at infinity or supersonic inflow, subsonic inflow associated with a jet 
or a nozzle, subsonic pressure outflow of given static pressure and supersonic 
outflow. Even if we still adopt a classical denomination for these boundary con- 
ditions, the fact that the fluid enters into the domain ("inflow") or goes outside 
it ("outflow") has no influence for the classification of the boundary conditions. 
The construction of a link between physical data and mathematical model is 
done via a precise choice of a boundary manifold. 

• Given state. This case corresponds typically to external aerodynamics 
problems. For this kind of flow, a state W r = Woo is known at a sufficient large 
distance between the object of study and the fluid boundary. In our present 
model this particular boundary is located at x = L. It is natural to set in a 
weak sense the boundary condition thanks to a (classical) Riemann problem and 
we have in this particular case a manifold M.^ te reduced to a single state : 

(3.4.1) M*£ e = {Woo}. 

We remark that this particular case is mathematically equivalent to the classical 
so-called "supersonic inflow" boundary condition. The knowledge of the entire 
state Woo is given and it interacts through the boundary with all the waves of 
the Riemann problem. 
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Figure 3.2 Subsonic jet inflow at the boundary of the domain {x < L}. 
It is described by a manifold .M ] q T parameterized by mass flux Q < and 
temperature T. 

• Subsonic jet inflow. In this case introduced in [CDV92] and described 
on Figure 3.2, we suppose that the fluid has a given flow rate Q and a given 
temperature T. We note that in order to respect the positive sign of the external 
normal of the boundary, the mass flux must be chosen negative (Q < 0). We 

introduce the boundary manifold M 3 q T equal to the set of states W that 
respect exactly the boundary condition : 

i<2 2 



(3.4.2) 



{ 



W= (p, q, e) 1 e Q, q = Q, e = pC v T + 



P 



} 



The resolution of the partial Riemann problem P(Wi, A4q t ) can be conducted 
in the plane of velocity and pressure as it is the case for the classical Riemann 
problem. We first express the internal energy e in terms of temperature T : 

(3.4.3) e = C V T 

and extract density p from mass flow Q and velocity u : 

(3.4.4) p=9.. 

u 

Then the thermostatic law (1.1.2) for polytropic perfect gas can be written in 
this context : 

(3.4.5) pu = (y-l)C v QT 

and the boundary manifold J\4q t is represented in velocity-pressure plane with 
an hyperbola, as shown on Figure 3.3. The resolution of the partial Riemann 
problem in space-time is easy : manifold J\4q t is of codimension 1 and from 

a given state Wi in the vicinity of J\Aq t , there exists a unique state W r 6 

■M 3 q T i issued from Wi through a 1-wave. In the case proposed on Figure 3.4, 
this wave is a shock wave and the effect of the interaction at the boundary is 
the incoming of a shock wave inside the domain {x < L}. 
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Figure 3.3 Resolution in velocity-pressure plane of the subsonic jet inflow 
boundary condition associated with the manifold M.q T - 




Figure 3.4 Resolution in space-time of the subsonic jet inflow boundary con- 
dition associated with the manifold M. ] q T - 

• Subsonic nozzle inflow. This case described on Figure 3.5. A given 

tank such the chamber of combustion of an engine to fix the ideas, contains some 
stagnation gas at rest at high temperature T s and high pressure p s . This gas is 
conducted to the inflow boundary through a compression wave in a convergent 
nozzle that maintains both total enthalpy H and specific entropy S. Recall 
that in a general manner, total enthalpy is defined according to 
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and specific entropy has been used at relation (1.1.12). When stagnation tem- 
perature T s and stagnation pressure p s are fixed, the total enthalpy H and 
the specific entropy £ are given by 

(3.4.7) H = C V T S1 E = C v Log( (^Y (■ 

\\T J \p J ) 

where Tq is the temperature of reference state in relation (1.1.12). We define 
the associated boundary manifold .M 1 ^ ff e as the set of states such that total 
enthalpy and entropy are equal respectively to H and £ ; we set 

( W = (p, q, e) 1 G Q, q = pu, 
(3 4 8) M™™ le = < e=\pu 2 + pe, p=(j-l)pe, 

„.2 , 



—u 



IP TT T ( PPQ \ S 

+ 7 TV" = H , L °g( r = 77" 




L, x 

Figure 3.5 Subsonic nozzle inflow at the boundary of the domain x2L. It is 
described by a manifold M.™ z ^ e parameterized by total enthalpy H and spe- 
cific entropy E and in an equivalent way by stagnation temperature T s and 
stagnation pressure p s . 

The "projection" of manifold .M£fff e in velocity-pressure plane is obtained by 
elimination of density p inside the relations presented at definition (3.4.8). It 
comes 

(3.4.9) -u 2 + -J—^l^/^p^ = H 

2 7-I/00 
and in terms of stagnation parameters : 
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The corresponding curve is represented on Figure 3.6 with the particular value 
of 7 = 



i.e. 



7" 



7 



2 
7' 



We have also represented a graphical resolution of the 
partial Riemann problem P(Wi, A^y^lf 16 ) in the same plane. 



CD shock 



nozzle 




max 



max 



Figure 3.6 Resolution in velocity-pressure plane of the subsonic nozzle inflow 
boundary condition associated with the manifold /V^|? le . 



• Subsonic pressure outflow. At a subsonic exit, it is classical to pre- 

scribe the static pressure. The linearized analysis shows that only one character- 
istic enters inside the domain of study and in consequence only one parameter 
is a priori necessary to close the boundary problem. The manifold associated to 
the pressure datum II is denoted by _/V4P ressure anc [ [ s particularity simple to 
define : 



(3.4.11) 



pressure 



W = (p,q = pu,e = pe + -p« 2 )^ G f2, 



(7 - l)pe 



n 



The resolution of the partial Riemann problem P(Wi, ^P ressure ) is presented 
on Figure 3.7. Note here that an inflow of the fluid is absolutly compatible with 
the manifold /Vt£[ ressure . 
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Figure 3.7 Subsonic pressure outflow at the boundary of the domain x2L. It 
is described by a manifold A^n essure . 

• Supersonic outflow. This case corresponds physically to natural boun- 

dary condition at the exit of a De Laval nozzle with sonic neck. No numerical 
datum has to be prescribed as we have seen previously with the linearized anal- 
ysis. Nevertheless, the fact that the limit state corresponds to a supersonic 
outflow (i.e. u > c) can enter in conflict with the initial datum for example. In 
[Du87] and [Du88], we have analyzed numerically this kind of problem. A nat- 
ural way to prescibe a "supersonic outflow" boundary condition is to introduce 
a boundary manifold ,M super that is now an half space as seen in section 2.5 : 

f TT7- / pc 2 1 2 \t ^ 

(3.4.12) M supei = I W = [P,Q = Pu,e= + ^pu ) eft, 

I u - c > 0. 

and to take into account the boundary condition with a partial Riemann problem 

P{Wi, .M super ) as seen previously at Proposition 2. 



3.5 Rigid wall and moving solid boundary. 

• The first important case in practice is the surface of some rigid body. In 
the multidimensional case, if n is the external normal to the body and u = 
(w, v) the bidimensional field of velocity, this boundary condition takes the form 
of an impermeability condition : 
(3.5.1) u.n = 0. 

At one space dimension, the external normal is equal to 1 and the condition 

(3.5.1) reduces to a nonlinear version of condition (3.2.11), id est 

(3.5.2) u = 0. 
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As in previous cases, this physically given boundary condition can be incom- 
patible with the state W(L~, t) present near the boundary. We propose here 
to write in a weak way the boundary condition (3.5.2) with a partial Riemann 

problem associated with boundary manifold M. r = .Mo el ° clty : 

(3.5.3) M? locity = {W=(p,q = pu,e) t GO, w = o}. 
Proposition 3. Stationary state for a rigid wall. 

Let W| be some given state, A^o el ° clty the manifold defined in (3.5.3) and 
W(x, t) = U(j ; Wi, 7Wo elocity ) the entropy solution of the partial Riemann 
problem P(W[, A^o el ° clty ) proposed at Theorem 2. Then the velocity u* of the 
stationary state U(0; Wi, .Mo elocity ) is nun an d we have 

(3.5.4) U(0 ; W h M v elocity ) = (p* , , J^f . 



(T)rarefaction 




U left <0 L 



right 



M 



(T) shock 
o 

velocity 



velocity 




u irft >0 L 



X 



Figure 3.8 Rigid wall boundary condition. The velocity of the stationary state 
is equal to zero. 

• The proof of Proposition 3 is clear if velocity ui of left state Wi is 

negative. In this case, the 1-wave that links state W\ and manifold A^o el ° clty 
is a rarefaction and the fan ]fi^ , fi^[ in space-time plane is defined by the 
conditions 

(3.5.5) ui — ci = jii < £ < (if = u r — c r = —c r < 
as presented on Figure 3.8. Then we have 

(3.5.6) U(0;W h Ml elocity ) = W r 

and property (3.5.4) is true in this particular case. If u\ > 0, the 1-wave is 
a shock, we have from (1.6.14) and (1.4.16) and due to the classical inequality 
p* > pi : 
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(3.5.7) 



o\ = ui - 




p* + fi 2 pi 



< ui - 



pi (p* + II 2 Pi) 



(p-pl) 



u 



Property (3.5.4) is a consequence of the general order (3.5.7) for the different 
waves : the celerity o\ of the 1-wave is less or equal than the celerity u* of the 
2- wave. For the partial Riemann problem P(Wi, .Mo elocity ), we have u* = 
and Proposition 3 is established. □ 

• A consequence of Proposition 3 is the fact that the only wave of par- 

tial Riemann problem P(Wi, A^o el ° clty ) is located in the quarter of space 
{x < L, t > 0} as depicted on Figure 3.8. We remark also that in this quar- 

velocity 



ter of space, the solution of the partial Riemann problem P(Wi, Ai 



) 



is 



identical to the solution of the Riemann problem R(Wi, Wi) between left state 

Wi and its "mirror" Wi. Recall that, following e.g. Roache [Ro72], the mirror 
state is defined by the conditions 

(3.5.8) pi = pi , ui = -ui , & = pi . 

One advantage of the notion of partial Riemamm problem is that we do not need 
to use this very numerical notion in what follows. 



(T) shock 



velocity 




w left 

(T) rarefaction 



Figure 3.9 Moving boundary condition at the boundary of the domain {x < 
L}. It is described by a manifold A4^ locity . 
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Figure 3.10 Solution in space-time of a moving boundary condition at the 
boundary with imposed velocity V. 

• The next example is motivated by fluid-structure interaction problems. 

In this kind of configuration, the boundary of the fluid domain is a moving 
boundary and the displacements are sufficiently small compared to typical 
length scales of the problem and can be neglected. Nevertheless, the boundary 
is moving and the value V of velocity is supposed to be given. We consider the 
weak treatment of the boundary condition u = V by the introduction of the 

manifold A ^velocity 



^veiocixy ^ e ^ ne( ^ as a generalization of (3.5.3) 



(3.5.9) 



^velocity 



U 



V 



} 



The resolution of the partial Riemann problem P(Wi, Aiy ° C1 y ) at the boun- 
dary is straightforward (see [Du99]). The unknown pressure p*(V) on the boun- 
dary is solution of the following equations introduced in (1.6.6) : 



f V - ui + V(p*(V); pi 

\ v - Ul + <p(p*(v) ]Pl 

that can be explicited in this case : 



(3.5.10) 



Pi 
Pi 



(3.5.11) 



V - ui + 



V - ^ + 



2c f 
7-1 



P*(V) 
Pi 



2 7 



7) = 
7) = 

- 1 












V2(p*(V)- Pl ) 



p*{V) < pi 
P*(V) > pi, 

if V - ui > 
= if V - ui < . 



f pi[(7 + l)p*(V) + (y-l)pi] 

With this approach, it is very easy to precise the difference p* (V) — p* (0) of 
boundary pressures associated with moving boundary with velocity V and a 
rigid wall. At the first order, we find : 

(3.5.12) p*(V) = p*(0) - piciV + 0{V 2 ). 

We have developed in [Du99] a complete procedure to take into account numeri- 
cally a displacement of the boundary of small amplitude and small celerity with 
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the so-called limiting flux for moving boundary. 



4) Application to the finite volume method. 
4.1 Godunov finite volume method. 

• The system of conservation laws 

(4.1.1) ftW(x, t) + -F(% t)) = 

is now discretized with the finite volume method. We consider a finite one- 
dimensional domain of study ]0, L[ and we divide it with a mesh T initially 
composed by a family of J+l vertices (note that integer J depends on 

the mesh T) : 

(4.1.2) S = { = xq < xi < ■ ■ ■ < Xj < Xj+i < • ■ ■ < xj = L } . 

We construct from vertices (%j)o<j<J the family £ of finite elements or 
control volume cells K : 

(4.1.3) ]0,L[= |J K, K=] Xj ,x j+1 [ e£ r , 0<j<J-l, 

KE£ T 

we remark that two different cells K\ and K2 have an intersection with null 
measure : 

(4.1.4) mes(i^ini^ 2 ) = 0, K 1 ^K 2 e£ r 

and we introduce the two vertices S-{K) and S+(K) that define the boundary 
dK of element K : 

(4.1.5) K — ]S-(K) , S + (K)[, OK = { S_(K) , S + (K)} , K e £^_ . 

• For each element K of the mesh [K G we introduce the mean 
value Wk of the solution of conservation law (4.1.1) : 

(4.1.6) W K = [ W(x)dx, Ke£ 

I K I Jk T 
and we consider the ordinary differential equation satisfied by the functions 

[0, +oo[ 3 1 1 — > Wxit) G ^- We integrate the conservation law (4.1.1) in space 

in each control volume K. After integrating by parts the divergent term ^—^p- , 
it comes 

(4.1.7) \K\ ^ +F (w(S + {K),t))-F(w(S-{K),t)) =0, K G £^_ . 

The expression (4.1.7) shows that a discretization procedure can be achieved if 
we are able to define the flux F$ for each vertex S G in terms of the data, 

i.e. the mean values {Wk , K G £ } and of the boundary conditions. This 
numerical modelling characterizes the so-called finite volume method (see e.g. 
Godunov et al [GZIKP79], Patankar [Pa80], Harten, Lax and Van Leer [HLV83] 
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or Faille, Gallouet and Herbin [FGH91] among others) and we will precise some 
efficient choices in practice in what follows. 

• We restrict ourselves for a time to vertices S E S that are such that 

T 

two finite elements Ki(S) and K r (S) possess the vertex S in their boundary. 
It is the case when vertex S is internal to the domain ]0, L[ : 

(4.1.8) dKi(S) ndK r (S) = {S} , S E <S^_ not located on the boundary. 

In practice, S is one of the vertices xi,-", xj_i and Ki(xj) =]xj-i, Xj[, 
K r (xj) =]xj , Xj+i[. With Godunov [Go59] we propose a numerical model for 
the numerical flux Fs at vertex S with the help of the Riemann problem. We 
denote by U(£ = f ; W\ , W r ) the entropy solution of the Riemann problem 
R(Wi, W r ). We set for internal vertices 

(4.1.9) F s = F(u(0-W Kl{s): W Kr{s) )y SE {*!,•••, 

and by this way, we have constructed an ordinary differential equation 

(4 - L10) + w\ ( Fs ^ K > ~ Fs -< K > ) = °- Ke£ r 

for control volumes K such that their boundary dK is not located on the 
boundary. 

4.2 Boundary fluxes. 

l W i r W. ' 

M left I 1 1 1 1 1 M right 

'x =0 x, ... x. , x. ... x T . x=L' 
I 1 j-i J J-1 J i 



Figure 4.1 Flux boundary conditions for the finite volume method. At the two 
boundaries of this unidimensional domain, the partial Riemann problems have 
to be set to compute the flux at the boundary. 

• We wish now to give a sense to equation (4.1.10) for all control volumes 
K of the mesh T . We make the hypothesis that at the left interface {x = 0} 
and at the right interface {x = L] of the domain, a left boundary manifold 
M-i and a right boundary manifold M. r are determined by the physical data 
(see Figure 4.1). Each boundary manifold is supposed to be chosen inside the 
family proposed previously at relations (3.4.1) for given state Woo, (3.4.2) for 
subsonic jet inflow M. ] q T (with mass flux parameter Q and temperature T), 
(3.4.8) for subsonic nozzle inflow .M^fff 6 (with numerical datum defined with 
total enthalpy H and specific entropy E), (3.4.11) for subsonic pressure outflow 
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^pressure ( w jth given static pressure II), (3.4.12) for supersonic outflow .A4 super 
or (3.5.9) for a moving boundary condition M™ locity (with celerity V) : 

(4.2.1) M i: M r e {{Woo} , M%\ T , M%g* , M% essure , .M super , A^ elocity } . 

• We denote by [/(£ = f ; Mi , W r ) and U(£ = f ; Wi , M r ) the self- 
similar entropy solution of the partial Riemann problems P(Mi , W r ) and 
P(Wi , M r ) proposed at Theorem 2 and relation (2.4.14). The boundary flux 
F XQ at the left interface is evaluated with the resolution of the partial Riemann 
problem P(Mi , W} XQjXl [) between the left boundary datum Mi and the state 
value W r in the first cell ]xq, x\\\ 

(4.2.2) F X0 = F(V(0;M, W ]xQtXl[ )) . 

In a similar way, the boundary flux F XJ at the right interface {x = L} is ob- 
tained thanks to the resolution of the partial Riemann problem P(W} Xj _ 1jXj [ , M 1 
between the left state W\ in the last cell ]xj-\, xj[ and the right boundary 
datum M. r '■ 

(4.2.3) F XJ = F(u(0-W ]xM[l M r )). 

With this choice of boundary fluxes F Xo and F XJ , the ordinary differential 
equation (4.1.10) takes a mathematical sense for all the control volumes K G £ . 

4.3 Strong nonlinearity at the boundary. 

• We focus here on the fact that the introduction of a partial Riemann pro- 
blem allows the treatment of strongly nonlinear effects at the boundary. Consider 
to fix the ideas the boundary {x = L} associated with computational domain 
]0, L[. We suppose that the physical conditions at this boundary are taken into 
account with the help of some manifold M. r . In order to consider weakly the 
boundary condition, we have introduced in (3.3.5) the set f3{M. r ) of admissible 
values at the boundary : 

o -i\ R( KA \ — / vames a t f = 0~ °f t ne entropic solution of the] 
[ } Py JVlr ) ~ \ partial Riemann problem P(W, M r ), W e £1 J ' 
Then the boundary condition at x = L is considered in the continuous case in 

(3.3.4) and we have by definition : 
(4.3.2) W(L~, t) e /3 (Mr) ■ 

This definition is compatible with the proposed implementation with the finite 
volume method : the boundary flux evaluated in (4.2.3) is the flux of the par- 
ticular state at £ = in the selfsimilar solution £/(£; Wi, M r ) of the partial 
Riemann problem P(Wi, M r ). 

• In lot of cases, the boundary condition acts as in the linear regime, 
i.e. the left state Wi in the partial Riemann problem P(Wi, M r ) lies in the 



53 



Francois Dubois 



vicinity of the manifold Ai r and the stationary value U(£ = 0; Wi, M. r ) is 
located inside the boundary manifold A4 r : 

(4.3.3) U(£ = 0; Wi, M. r ) £ M-r : weak nonlinearity at the boundary. 
This case occurs typically for a manifold of codimension 1. Consider to fix the 
ideas a manifold _/v^P ressure where the static pressure is fixed by the physical 
outflow conditions. The linearized approach in the vicinity of the (unknown) 
state W* = U(£ — 0; Wi, M. r ) is in general correct (see Figure 4.2) and 
the classical approach (see e.g. Viviand and Veuillot [VV78], Chakravarthy 
[Ch83], Osher and Chakravarthy [OC83]) consists in determining a state W* 
at the boundary that satifies strongly the boundary condition, i.e. satisfies the 
condition W* G M r and satisfies also the linearized gas dynamics equations 
written under a characteristic form along the two characteristics going outside 
the domain of study. 



A t 

u* 



u* — c* 



1 -wave 




W 



1 2-wave 



u*+ c* 



3 -wave 



/ ^ 



x 



X 



Figure 4.2 Weakly nonlinear boundary condition. The state W* at the boun- 
dary x=x* belongs to the boundary manifold. 

• This approach is in defect when the limit state is far from the boun- 

dary manifold M. r . By example, a supersonic boundary condition is prescribed 
but the limiting state W(L~,t) corresponds to a subsonic outflow. On the 
contrary, a subsonic outflow with pressure II is supposed to be given but the 
limiting state W(L~ : t) is associated with a supersonic exit ! The flexibility of 
the partial Riemann problem allows the treatment with strong nonlinear waves. 
When condition (4.3.3) is in defect, the boundary condition (4.2.2) or (4.2.3) 
takes into account the physical data that constructs the interaction at the boun- 
dary whereas these data are not directly used for the final computation of the 
boundary flux (4.2.2) or (4.2.3). 
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4.4 Extension to second order accuracy and to two space dimensions. 

• We detail in this section a generalization for unstructured meshes of the 
"Multidimensional Upwindcentered Scheme for Conservation Laws" proposed by 
Van Leer [VL79]. At one space dimension on a uniform mesh, it is classical to 
consider a scalar field z among the primitive variables, i.e. 

(4.4.1) z E {p,u,v,p} (primitive variables) 

and instead of computing the interface flux with relation (4.1.9), to first construct 
two interface states W<7 and W$ on each side of the interface S. Then the 
flux is evaluated by the decomposition of the discontinuity : 

(4.4.2) F S = F(u{0;Ws,W+)), S E { x u ■ ■ ■ , xj- X } . 

This nonlinear interpolation is done with a so-called "slope limiter" ip(») that 
operates on each variable proposed in (4.4.1) and we have typically when a 
left-right invariance is assumed [Du91] : 



(4.4.3) 
(4.4.4) 




Figure 4.3 Structured cartesian mesh. The control volumes are exactly the 
elements of mesh T. 
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Figure 4.4 Unstructured mesh composed by triangular elements. The control 
volumes are exactly the elements of mesh T. 



• We focus now on two problems for the extension to second order ac- 

curacy of Godunov finite volume method. One is associated with the use of 
unstructured meshes and the other with the treatment of boundary conditions. 
As in the one-dimensional case, the domain of study is decomposed into finite 
elements (or control volumes) K G 8 than can be structured in a cartesian 

way (Figure 4.3) or with a cellular complex as in Figure 4.4. In both cases, the 
intersection of two finite elements define an interface / G . We denote by 

iif the normal at the interface / that separates a left control volume Ki(f) 
and a right control volume K r (f). The ordinary differential equation (4.1.7) is 
replaced by a multidimensional version : 

(4-4.5) \K\ ^ + J2 |/| <$>(W K , n f , W Kr{f) ) =0, K G ^ . 

fCdK 

For internal interfaces, the function , n/ , •) is equal to the flux of the so- 
lution at j = of the Riemann problem between states W^m and Wx r (f) 
in the one-dimensional direction along normal n/ in order to take into ac- 
count the invariance by rotation of the equations of gas dynamics (see [GR96]). 
We suppose also that for each interface / of the boundary, a boundary man- 
ifold A4f of codimension p(f) is given and the normal direction iif is by 
convention external to the domain of study. In consequence, the control cell 
K = Ki(f) has the face / in its boundary (/ C d(Ki(f))) and in relation 
(4.4.5), the state Wx r {f) belongs to A4f {Wx r (f) EM.f) and is equal to the 
state W p(y ^ introduced in relation (2.4.2) when solving the partial Riemann 
problem P(Wxj(/) , n/ , M.f) between state Wxdf) and manifold M.f- 
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Figure 4.5 Cellular complex mesh with triangles and quadrangles. Three 
neighbouring cells are necessary to determine the gradient in triangle K and to 
limit eventually its variation. 

• We consider now a finite element K internal to the domain. The exten- 
sion to second order accuracy of the finite volume scheme consists in replacing 
the arguments W Kl ^ and Wx r (f) in relation (4.4.5) by nonlinear extrapola- 
tions Wx t (f),f and Wx r (f),f on each side of the boundary of state data and 
evaluated as described in what follows. We first introduce the set J\f(K) of 
neighbouring cells of given finite element K E 8 , as illustrated on Figure 4.5 : 

(4.4.6) Af(K) = {LeS^ ]/G7 r , / C OK n dL } . 

For L E J\f(K), we suppose by convention that the normal rif to the face 
/ C dK n dL is external to the element K id est K r {f) = K, Ki(f) = L. 
We introduce also the point VK,f on the interface / C dK that links the 
barycenters xk and x Kr ^ : 

(4_4_ 7 ) / VK,f = {l-0 K ,f)x K + 6 K , f x Kr(f) , y K ,f£f, 

y f C dK , K finite element internal to mesh T . 

Then, following Pollet [Po88], for z equal to one scalar variable of the family : 

(4.4.8) z E {p,pu,pv,p} 

we evaluate a mean value zk, J on the interface / : 

(4.4.9) zwj = {i-e K , f )z K + e KJ z KrU) 

and the gradient Vz(K) of field z(») in volume K with a Green formula : 

* ndl= TK~\ ^ \f\z]^Jn f ,KeS T . 

1 1 fcdK 



(4.4.10) Vz{K) = 
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• An ideal extrapolation of field z(») at the interface / would be : 

(4.4.11) z K ,f = z K + Vz(K).(y K ,f - x K ) 

but the corresponding scheme is unstable as seen by Van Leer [VL77]. When 
the variation Vz(K) • (yx,f — %k) is too important, it has to be "limited" as 
first suggested by Van Leer [VL77]. For doing this in a very general way, we 
introduce the minimum 171k{z) and the maximum Mr(z) of field z(») in the 
neighbouring cells : 

(4.4.12) m K (z) = min{z L , LeAf(K)} 

(4.4.13) M K {z) = max{z L , L G M{K) } . 
If the value zk is extremum among the neighbouring ones, i.e. if 
or zk > Mk(z), we impose that the interpolated value zx,f 
cell value zk '• 

(4.4.14) z K ,f = z K if z K < m K {z) or z K > M K {z) , / C dK . 
When on the contrary zk lies inside the interval [rriK(z), Mk(z)], we impose 
that the variation zk, / — zk is limited by some coefficient k (0 < k < 1) 
of the variations zk — ttik{z) and Mk(z) — zk- We introduce a nonlinear 
extrapolation of the field z(») between center xk and boundary face yx / (/ C 
dK) : 

(4.4.15) z K , f = z K + a K {z)Vz{K).(y KJ -x K ) , f C dK 
with a limiting coefficient olk{z) satisfying the following conditions : 

< olk{z) < 1 , z{») scalar field defined in (4.4.8) 



zk < rn K {z), 
is equal to the 



K G S 



T 



(4.4.16) <^ k(z K -m K (z)) < a K (z)Vz(K) . (y Ki f - x K ) < k{M K {z)-z K ) 



VfcdK, Ke£. 



T 




Figure 4.6 Examples of limiter functions that can be easily extended to un- 
structured meshes. 

Then olk{z) is chosen as big as possible and inferior or equal to 1 in order to 
satisfy the constraints (4.4.16) as displayed on Figure 4.6 : 
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(AA ^s i \ • r, 7 mm(M K (z) -z K , z K -m K {z)) 

(4.4.17) a K (z) = mm l,k TrvTT^\ < w t 7- 

max { | Vz(K) • {y K) f - x K ) \ , / C dK ) 

• In the one dimensional case with a regular mesh, it is an exercice to 

re- write the extrapolation (4.4.15) under the usual form (4.4.3) in the context 
of finite differences. In this particular case, some limiter functions r i — > ipk(r) 
associated with particular parameters k are shown on Figure 4.6. For k — 1, 
we recover the initial limiter proposed by Van Leer in the fourth paper of the 
family "Towards the ultimate finite difference scheme..." [VL77] ; for this reason, 
we have named it the "Towards 4" limiter (see Figure 4.6). When k — \ we 
obtain the "min-mod" limiter proposed by Harten [Ha83]. The intermediate 
value k = | is a good compromise between the "nearly unstable" choice k — 1 
and the "too compressive" min-mod choice. We have named it STS and it has 
been chosen for our Euler computations in [DM92]. 




Figure 4.7 Slope limitation at a fluid boundary. 

• We explain now the way the preceeding scheme is adapted near the 

boundary. We first consider a fluid boundary. When K is a finite element 
with some face g C dK lying on the boundary, we still define the set J\f(K) of 
neighbouring cells by the relation (4.4.6) as shown on Figure 4.7. The number 
of neighbouring cells is just less important in this case. Then points UK,f 
are introduced by relation (4.4.7) if face / does not lie on the boundary and 
by taking the barycenter of face g if it is lying on the boundary. The only 
difference is the way the values zk, g are extrapolated for the face g that is 
on the boundary ; we set 

(4.4.18) zk, g = zk , g C dK , g face lying on the boundary of the do- 
main. 
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When values zk, / are determined for all the faces / C dK, the gradient 
V 'z(K), the minimal ttik^z) and maximal Mk{z) values among the neigh- 
bouring cells are still determined with the relations (4.4.10), (4.4.12) and (4.4.13) 
respectively. The constraints (4.4.16) remain unchanged except that no limita- 
tion process is due to the faces lying on the boundary. In a precise way, we 
set : 

k min( Mk{z) — zk , %k ~ m K{z) ) 
' max { | VzjK) . (y K , f - x K ) | , / C dK, K r (f) eAf(K) } _ ' 
Then the interpolated values zk, / for all the faces / C dK are again predicted 
with the help of relation (4.4.15). 

• For a rigid wall, the limitation process is a little modified, as presented 
at Figure 4.8. First we introduce the limit face g inside the set of neighbouring 
cells : 



(4.4.20) 




Figure 4.8 Slope limitation at a solid boundary. 

For the face(s) g C dK lying on the solid boundary, we determine preliminary 
values zk, g by taking in consideration at this level the nonpenetrability boun- 
dary condition (3.5.1). We introduce the two components n x g and n v g of the 
normal n g at the boundary and we set, in coherence with variables (4.4.8) : 



(4.4.19) a K {z)=mm 
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PK,g = PK 

(4 4 21) { uW^ = p K {uK-(uK»n g )n x g ) 

| ~pKr g 7 WTg = Pk {vk ~ {u K •n g )n y g ) 
V~K~g = PK ■ 

We consider also these values for the limitation algorithm. We define "external 
values" zl for L = g and face g lying on the boundary as equal to the ones 
defined in relation (4.4.21) : 

(4.4.22) z g = ^K~gi z {*) field defined in (4.4.21), g C dK on the boundary. 
Then the extrapolation algorithm that conducts to relation (4.4.15) for extrap- 
olated values zk, / is used as in the internal case. 

• When all values zk, / are known for all control volumes K £ £ , all 

faces / CK and all fields z(») defined at relation (4.4.8), extrapolated states 
Wk, f are naturally defined by going back to the conservative variables. Then 
we introduce these states as arguments of the flux function $(•, n/ , •) and 
obtain by this way a new system of ordinary differential equations : 

(4.4.23) \K\ ^ + E l/l HW K ,f,n f ,W Krif) , f ) =0, K £ £^_ . 

fCdK 

The numerical integration of such kind of system is just a question of Runge- 
Kutta scheme as presented in [CDV92]. We have used with success in [DM92] 
the Heun scheme of second order accuracy for discrete integration of (4.4.23) 
between time steps n At and (n + 1) At : 



(4.4.24) Li (W K -WZ)+Y, l/l ^{ w lf^f^ W K A f),f) Ke£ 

fcdK 

(4.4.25) \^[Wk-Wk)+Y. l/l *(w K ,f,n f ,W Kr{f) , f ) = 0, K £ ^ 



fCdK 



(4.4.26) W£ +1 = - (W K + W%) , K £ ^ . 
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